#Transforming roll-call data & modelling
rm(list = ls()) #clear environment
library(tidyverse)
library(glue)
library(data.table)
library(DataCombine)
library(MCMCpack)
library(coda)
library(car)
library(dplyr)
library(rjags)
library(extrafont)
p <- list.functions.in.file("IRT_models.R")
p
#Transforming roll-call data & modelling
rm(list = ls()) #clear environment
library(tidyverse)
library(coda)
library(car)
library(dplyr)
library(ggplot2)
library(rjags)
######################################
#Import and organise divisions data
######################################
#Import MP names
mpnames <- read.delim("./Other data sources/mps-2017.txt", header = TRUE, sep = "\t")
mpnames <- unite(mpnames, mpname, c(firstname, surname), remove=TRUE, sep = " ") #merge first name and surname
mpnames <- mpnames[,c(1,2)] #keep only MP name and id
mpnames$mpname <- gsub("O&#39;", "O'", mpnames$mpname) #fix encoding issue
#Import roll call matrix
#transform into data frame with votenumbers as row names, mp names as column names
rollmatrix <- read.delim("./Other data sources/votematrix-2017.txt", header = TRUE, sep = "\t")
rollmatrix$Bill <- gsub("&#8217;", "'", rollmatrix$Bill) #fix encoding issue
rollmatrix$date <- as.character(rollmatrix$date)
rollmatrix$voteno <- as.character(rollmatrix$voteno)
#Select Brexit-related bills
#see document called BrexitBillSelect.docx to understand which bills are deemed relevant.
inxvoteno <- grep("293|307|308|310|311|312|331|332|333|345|346|347|354|357|358|359|360|361|362|363|364|374|386|387|388|389|390|391|392|393|394|395|397|398|399|400|404|405|406|407|408|409|439|440|441|442|\\<3\\>|\\<4\\>|\\<5\\>|\\<11\\>", rollmatrix$voteno, fixed = F)
rollmatrix <- rollmatrix[inxvoteno,]
rollmatrix <- rollmatrix[rollmatrix$date > as.Date("2019-01-01"),]
rownames(rollmatrix) <- rollmatrix$voteno # realign row names
rolls <- rollmatrix[,-c(1,2,3,4)] #drop some columns (rowid, date, voteno, bill)
rolls <- subset(rolls, select=-X) #remove this empty column named 'X'
#Transform vote matrix
colnames(rolls) <- gsub("mpid", "", colnames(rolls)) #remove "mpid" from variables
rolls <- t(rolls) #transpose matrix
#Replace voting codes
rolls <- recode(rolls, "2=10") #yes
rolls <- recode(rolls, "1=10") #yes
rolls <- recode(rolls, "-9=0") #NA
rolls <- recode(rolls, "3=0") #NA
rolls <- recode(rolls, "4=20") #no
rolls <- recode(rolls, "5=20") #no
#Join MPs' names
rolls <- cbind(rolls, mpid = unlist(rownames(rolls)))
rolls <- merge(mpnames, rolls, by = "mpid", all.y = TRUE) #match mp names with mp ids
rolls <- cbind(rolls[,1:2], mutate_all(rolls[,3:52], function(x) if(is.character(x)) as.numeric(as.character(x)) else x)) #change votes back to numeric
#Collapse data by MP (for MPs who switched party affiliation)
rolls <- rolls %>%
group_by(mpname) %>%
summarise_all(list(max)) #summarize voting records by mp (deprecated function warning but works just fine)
#Recode again to a 0-1-NA format necessary for IRT models
rolls <- replace(rolls, rolls == 0, NA) #NA
rolls <- replace(rolls, rolls == 10, 1) #YES
rolls <- replace(rolls, rolls == 20, 0) #NO
#Get rid of MPs with huge gaps in voting record
rolls$na_count_p1 <- apply(rolls[,11:52], 1, function(x) sum(is.na(x))) #no of missing votes in period 1
rolls$na_count_p2 <- apply(rolls[,3:10], 1, function(x) sum(is.na(x))) #no of missing votes in period 2
rolls <- rolls[rolls$na_count_p1 > 21 | rolls$na_count_p2 > 4,] #exclude if more than half of votes in either period missing (applies to 13 MPs)
#Identify anchors
rolls$anchorMP <- NA
rolls$anchorMP[rolls$mpname == "Steven Baker"] <- 1
rolls$anchorMP[rolls$mpname == "David Lammy"] <- 2
rolls <- rolls[order(rolls$anchorMP),-c(53)] #put anchor MPs to the top
#Restrict to Ys
rolls.jags <- rolls[,3:52]
rolls.jags <- as.matrix(rolls.jags)
###########################################
#Dynamic IRT Model with Quadratic Effects
###########################################
#JAGS code for model with fixed evolution parameter
dyn_irt_model_quadratic <- 'model{
#loop through actors
for(i in 1:nactors){
#loop through votes
for(j in 1:nvotes){
Y[i,j] ~ dbern(mu[i,j])
logit(mu[i,j]) <- alpha[j] + beta[j] * theta[actor[i], period[j]] + beta2[j] * theta[actor[i], period[j]] * theta[actor[i], period[j]]
}
change[i] <- theta[i,2] - theta[i,1]
}
#set normal priors on betas and alphas
for(j in 1:nvotes){
beta[j] ~ dnorm(0, 1)
beta2[j] ~ dnorm(0, 1)
alpha[j] ~ dnorm(0, 1)
}
#set priors on thetas (fix space with two actors)
theta[1, 1] ~ dnorm(1, 1)
theta[2, 1] ~ dnorm(-1, 1)
for(t in 2:nperiods){
theta[1, t] ~ dnorm(theta[1, t-1], tau.evol)
theta[2, t] ~ dnorm(theta[2, t-1], tau.evol)
}
for(c in 3:nactors){
theta[c, 1] ~ dnorm(0, 1)
for(t in 2:nperiods){
theta[c, t] ~ dnorm(theta[c, t-1], tau.evol)
}
}
#set priors on evolution parameters for thetas
tau.evol <- 1
}
'
#Set parameters
jags.data.irt <- list(
Y = rolls.jags,
nactors = 630,
nvotes = 50,
nperiods = 2,
period = c(2, 2, 2, 2, 2, 2, 2, 2, rep(1, 42)),
actor = as.numeric(c(1:630))
)
#Fit JAGS model
model.jags.irt <- jags.model(file=textConnection(dyn_irt_model_quadratic), data = jags.data.irt, n.chains = 1, inits=list(.RNG.name="base::Wichmann-Hill",.RNG.seed=2355328))
#Transforming roll-call data & modelling
rm(list = ls()) #clear environment
library(tidyverse)
library(glue)
library(data.table)
library(DataCombine)
library(MCMCpack)
library(coda)
library(car)
library(dplyr)
library(rjags)
library(extrafont)
######################################
#Import and organise divisions data
######################################
#Import MP names
mpnames <- read.delim("./Other data sources/mps-2017.txt", header = TRUE, sep = "\t")
mpnames <- unite(mpnames, mpname, c(firstname, surname), remove=TRUE, sep = " ") #merge first name and surname
mpnames <- mpnames[,c(1,2)] #keep only MP name and id
mpnames$mpname <- gsub("O&#39;", "O'", mpnames$mpname) #fix encoding issue
#Import roll call matrix
#transform into data frame with votenumbers as row names, mp names as column names
rollmatrix <- read.delim("./Other data sources/votematrix-2017.txt", header = TRUE, sep = "\t")
rollmatrix$Bill <- gsub("&#8217;", "'", rollmatrix$Bill) #fix encoding issue
rollmatrix$date <- as.character(rollmatrix$date)
rollmatrix$voteno <- as.character(rollmatrix$voteno)
#Select Brexit-related bills
#see document called BrexitBillSelect.docx to understand which bills are deemed relevant.
inxvoteno <- grep("293|307|308|310|311|312|331|332|333|345|346|347|354|357|358|359|360|361|362|363|364|374|386|387|388|389|390|391|392|393|394|395|397|398|399|400|404|405|406|407|408|409|439|440|441|442|\\<3\\>|\\<4\\>|\\<5\\>|\\<11\\>", rollmatrix$voteno, fixed = F)
rollmatrix <- rollmatrix[inxvoteno,]
rollmatrix <- rollmatrix[rollmatrix$date > as.Date("2019-01-01"),]
rownames(rollmatrix) <- rollmatrix$voteno # realign row names
rolls <- rollmatrix[,-c(1,2,3,4)] #drop some columns (rowid, date, voteno, bill)
rolls <- subset(rolls, select=-X) #remove this empty column named 'X'
#Transform vote matrix
colnames(rolls) <- gsub("mpid", "", colnames(rolls)) #remove "mpid" from variables
rolls <- t(rolls) #transpose matrix
#Replace voting codes
rolls <- recode(rolls, "2=10") #yes
rolls <- recode(rolls, "1=10") #yes
rolls <- recode(rolls, "-9=0") #NA
rolls <- recode(rolls, "3=0") #NA
rolls <- recode(rolls, "4=20") #no
rolls <- recode(rolls, "5=20") #no
#Join MPs' names
rolls <- cbind(rolls, mpid = unlist(rownames(rolls)))
rolls <- merge(mpnames, rolls, by = "mpid", all.y = TRUE) #match mp names with mp ids
rolls <- cbind(rolls[,1:2], mutate_all(rolls[,3:52], function(x) if(is.character(x)) as.numeric(as.character(x)) else x)) #change votes back to numeric
#Collapse data by MP (for MPs who switched party affiliation)
rolls <- rolls %>%
group_by(mpname) %>%
summarise_all(list(max)) #summarize voting records by mp (deprecated function warning but works just fine)
#Recode again to a 0-1-NA format necessary for IRT models
rolls <- replace(rolls, rolls == 0, NA) #NA
rolls <- replace(rolls, rolls == 10, 1) #YES
rolls <- replace(rolls, rolls == 20, 0) #NO
#Get rid of MPs with huge gaps in voting record
rolls$na_count_p1 <- apply(rolls[,11:52], 1, function(x) sum(is.na(x))) #no of missing votes in period 1
rolls$na_count_p2 <- apply(rolls[,3:10], 1, function(x) sum(is.na(x))) #no of missing votes in period 2
rolls <- rolls[rolls$na_count_p1 > 21 | rolls$na_count_p2 > 4,] #exclude if more than half of votes in either period missing (applies to 13 MPs)
#Identify anchors
rolls$anchorMP <- NA
rolls$anchorMP[rolls$mpname == "Steven Baker"] <- 1
rolls$anchorMP[rolls$mpname == "David Lammy"] <- 2
rolls <- rolls[order(rolls$anchorMP),-c(53)] #put anchor MPs to the top
#Restrict to Ys
rolls.jags <- rolls[,3:52]
rolls.jags <- as.matrix(rolls.jags)
###########################################
#Dynamic IRT Model with Quadratic Effects
###########################################
#JAGS code for model with fixed evolution parameter
dyn_irt_model_quadratic <- 'model{
#loop through actors
for(i in 1:nactors){
#loop through votes
for(j in 1:nvotes){
Y[i,j] ~ dbern(mu[i,j])
logit(mu[i,j]) <- alpha[j] + beta[j] * theta[actor[i], period[j]] + beta2[j] * theta[actor[i], period[j]] * theta[actor[i], period[j]]
}
change[i] <- theta[i,2] - theta[i,1]
}
#set normal priors on betas and alphas
for(j in 1:nvotes){
beta[j] ~ dnorm(0, 1)
beta2[j] ~ dnorm(0, 1)
alpha[j] ~ dnorm(0, 1)
}
#set priors on thetas (fix space with two actors)
theta[1, 1] ~ dnorm(1, 1)
theta[2, 1] ~ dnorm(-1, 1)
for(t in 2:nperiods){
theta[1, t] ~ dnorm(theta[1, t-1], tau.evol)
theta[2, t] ~ dnorm(theta[2, t-1], tau.evol)
}
for(c in 3:nactors){
theta[c, 1] ~ dnorm(0, 1)
for(t in 2:nperiods){
theta[c, t] ~ dnorm(theta[c, t-1], tau.evol)
}
}
#set priors on evolution parameters for thetas
tau.evol <- 1
}
'
#Set parameters
jags.data.irt <- list(
Y = rolls.jags,
nactors = 630,
nvotes = 50,
nperiods = 2,
period = c(2, 2, 2, 2, 2, 2, 2, 2, rep(1, 42)),
actor = as.numeric(c(1:630))
)
#Fit JAGS model
model.jags.irt <- jags.model(file=textConnection(dyn_irt_model_quadratic), data = jags.data.irt, n.chains = 1, inits=list(.RNG.name="base::Wichmann-Hill",.RNG.seed=2355328))
#Transforming roll-call data & modelling
rm(list = ls()) #clear environment
library(tidyverse)
library(glue)
library(data.table)
library(DataCombine)
library(MCMCpack)
library(coda)
library(car)
library(dplyr)
library(rjags)
library(extrafont)
######################################
#Import and organise divisions data
######################################
#Import MP names
mpnames <- read.delim("./Other data sources/mps-2017.txt", header = TRUE, sep = "\t")
mpnames <- unite(mpnames, mpname, c(firstname, surname), remove=TRUE, sep = " ") #merge first name and surname
mpnames <- mpnames[,c(1,2)] #keep only MP name and id
mpnames$mpname <- gsub("O&#39;", "O'", mpnames$mpname) #fix encoding issue
#Import roll call matrix
#transform into data frame with votenumbers as row names, mp names as column names
rollmatrix <- read.delim("./Other data sources/votematrix-2017.txt", header = TRUE, sep = "\t")
rollmatrix$Bill <- gsub("&#8217;", "'", rollmatrix$Bill) #fix encoding issue
rollmatrix$date <- as.character(rollmatrix$date)
rollmatrix$voteno <- as.character(rollmatrix$voteno)
#Select Brexit-related bills
#see document called BrexitBillSelect.docx to understand which bills are deemed relevant.
inxvoteno <- grep("293|307|308|310|311|312|331|332|333|345|346|347|354|357|358|359|360|361|362|363|364|374|386|387|388|389|390|391|392|393|394|395|397|398|399|400|404|405|406|407|408|409|439|440|441|442|\\<3\\>|\\<4\\>|\\<5\\>|\\<11\\>", rollmatrix$voteno, fixed = F)
rollmatrix <- rollmatrix[inxvoteno,]
rollmatrix <- rollmatrix[rollmatrix$date > as.Date("2019-01-01"),]
rownames(rollmatrix) <- rollmatrix$voteno # realign row names
rolls <- rollmatrix[,-c(1,2,3,4)] #drop some columns (rowid, date, voteno, bill)
rolls <- subset(rolls, select=-X) #remove this empty column named 'X'
#Transform vote matrix
colnames(rolls) <- gsub("mpid", "", colnames(rolls)) #remove "mpid" from variables
rolls <- t(rolls) #transpose matrix
#Replace voting codes
rolls <- recode(rolls, "2=10") #yes
rolls <- recode(rolls, "1=10") #yes
rolls <- recode(rolls, "-9=0") #NA
rolls <- recode(rolls, "3=0") #NA
rolls <- recode(rolls, "4=20") #no
rolls <- recode(rolls, "5=20") #no
#Join MPs' names
rolls <- cbind(rolls, mpid = unlist(rownames(rolls)))
rolls <- merge(mpnames, rolls, by = "mpid", all.y = TRUE) #match mp names with mp ids
rolls <- cbind(rolls[,1:2], mutate_all(rolls[,3:52], function(x) if(is.character(x)) as.numeric(as.character(x)) else x)) #change votes back to numeric
#Collapse data by MP (for MPs who switched party affiliation)
rolls <- rolls %>%
group_by(mpname) %>%
summarise_all(list(max)) #summarize voting records by mp (deprecated function warning but works just fine)
#Recode again to a 0-1-NA format necessary for IRT models
rolls <- replace(rolls, rolls == 0, NA) #NA
rolls <- replace(rolls, rolls == 10, 1) #YES
rolls <- replace(rolls, rolls == 20, 0) #NO
#Get rid of MPs with huge gaps in voting record
rolls$na_count_p1 <- apply(rolls[,11:52], 1, function(x) sum(is.na(x))) #no of missing votes in period 1
rolls$na_count_p2 <- apply(rolls[,3:10], 1, function(x) sum(is.na(x))) #no of missing votes in period 2
rolls <- rolls[rolls$na_count_p1 > 21 | rolls$na_count_p2 > 4,] #exclude if more than half of votes in either period missing (applies to 13 MPs)
#Identify anchors
rolls$anchorMP <- NA
rolls$anchorMP[rolls$mpname == "Steven Baker"] <- 1
rolls$anchorMP[rolls$mpname == "David Lammy"] <- 2
rolls <- rolls[order(rolls$anchorMP),-c(53)] #put anchor MPs to the top
#Restrict to Ys
rolls.jags <- rolls[,3:52]
rolls.jags <- as.matrix(rolls.jags)
###########################################
#Dynamic IRT Model with Quadratic Effects
###########################################
#JAGS code for model with fixed evolution parameter
dyn_irt_model_quadratic <- 'model{
#loop through actors
for(i in 1:nactors){
#loop through votes
for(j in 1:nvotes){
Y[i,j] ~ dbern(mu[i,j])
logit(mu[i,j]) <- alpha[j] + beta[j] * theta[actor[i], period[j]] + beta2[j] * theta[actor[i], period[j]] * theta[actor[i], period[j]]
}
change[i] <- theta[i,2] - theta[i,1]
}
#set normal priors on betas and alphas
for(j in 1:nvotes){
beta[j] ~ dnorm(0, 1)
beta2[j] ~ dnorm(0, 1)
alpha[j] ~ dnorm(0, 1)
}
#set priors on thetas (fix space with two actors)
theta[1, 1] ~ dnorm(1, 1)
theta[2, 1] ~ dnorm(-1, 1)
for(t in 2:nperiods){
theta[1, t] ~ dnorm(theta[1, t-1], tau.evol)
theta[2, t] ~ dnorm(theta[2, t-1], tau.evol)
}
for(c in 3:nactors){
theta[c, 1] ~ dnorm(0, 1)
for(t in 2:nperiods){
theta[c, t] ~ dnorm(theta[c, t-1], tau.evol)
}
}
#set priors on evolution parameters for thetas
tau.evol <- 1
}
'
View(rolls)
#Transforming roll-call data & modelling
rm(list = ls()) #clear environment
library(tidyverse)
library(coda)
library(car)
library(dplyr)
library(ggplot2)
library(rjags)
######################################
#Import and organise divisions data
######################################
#Import MP names
mpnames <- read.delim("./Other data sources/mps-2017.txt", header = TRUE, sep = "\t")
mpnames <- unite(mpnames, mpname, c(firstname, surname), remove=TRUE, sep = " ") #merge first name and surname
mpnames <- mpnames[,c(1,2)] #keep only MP name and id
mpnames$mpname <- gsub("O&#39;", "O'", mpnames$mpname) #fix encoding issue
#Import roll call matrix
#transform into data frame with votenumbers as row names, mp names as column names
rollmatrix <- read.delim("./Other data sources/votematrix-2017.txt", header = TRUE, sep = "\t")
rollmatrix$Bill <- gsub("&#8217;", "'", rollmatrix$Bill) #fix encoding issue
rollmatrix$date <- as.character(rollmatrix$date)
rollmatrix$voteno <- as.character(rollmatrix$voteno)
#Select Brexit-related bills
#see document called BrexitBillSelect.docx to understand which bills are deemed relevant.
inxvoteno <- grep("293|307|308|310|311|312|331|332|333|345|346|347|354|357|358|359|360|361|362|363|364|374|386|387|388|389|390|391|392|393|394|395|397|398|399|400|404|405|406|407|408|409|439|440|441|442|\\<3\\>|\\<4\\>|\\<5\\>|\\<11\\>", rollmatrix$voteno, fixed = F)
rollmatrix <- rollmatrix[inxvoteno,]
rollmatrix <- rollmatrix[rollmatrix$date > as.Date("2019-01-01"),]
rownames(rollmatrix) <- rollmatrix$voteno # realign row names
rolls <- rollmatrix[,-c(1,2,3,4)] #drop some columns (rowid, date, voteno, bill)
rolls <- subset(rolls, select=-X) #remove this empty column named 'X'
#Transform vote matrix
colnames(rolls) <- gsub("mpid", "", colnames(rolls)) #remove "mpid" from variables
rolls <- t(rolls) #transpose matrix
#Replace voting codes
rolls <- recode(rolls, "2=10") #yes
rolls <- recode(rolls, "1=10") #yes
rolls <- recode(rolls, "-9=0") #NA
rolls <- recode(rolls, "3=0") #NA
rolls <- recode(rolls, "4=20") #no
rolls <- recode(rolls, "5=20") #no
#Join MPs' names
rolls <- cbind(rolls, mpid = unlist(rownames(rolls)))
rolls <- merge(mpnames, rolls, by = "mpid", all.y = TRUE) #match mp names with mp ids
rolls <- cbind(rolls[,1:2], mutate_all(rolls[,3:52], function(x) if(is.character(x)) as.numeric(as.character(x)) else x)) #change votes back to numeric
#Collapse data by MP (for MPs who switched party affiliation)
rolls <- rolls %>%
group_by(mpname) %>%
summarise_all(list(max)) #summarize voting records by mp (deprecated function warning but works just fine)
#Recode again to a 0-1-NA format necessary for IRT models
rolls <- replace(rolls, rolls == 0, NA) #NA
rolls <- replace(rolls, rolls == 10, 1) #YES
rolls <- replace(rolls, rolls == 20, 0) #NO
#Get rid of MPs with huge gaps in voting record
rolls$na_count_p1 <- apply(rolls[,11:52], 1, function(x) sum(is.na(x))) #no of missing votes in period 1
rolls$na_count_p2 <- apply(rolls[,3:10], 1, function(x) sum(is.na(x))) #no of missing votes in period 2
rolls <- rolls[rolls$na_count_p1 <= 21 & rolls$na_count_p2 <= 4,] #exclude if more than half of votes in either period missing
#Identify anchors
rolls$anchorMP <- NA
rolls$anchorMP[rolls$mpname == "Steven Baker"] <- 1
rolls$anchorMP[rolls$mpname == "David Lammy"] <- 2
rolls <- rolls[order(rolls$anchorMP),-c(53)] #put anchor MPs to the top
#Restrict to Ys
rolls.jags <- rolls[,3:52]
rolls.jags <- as.matrix(rolls.jags)
###########################################
#Dynamic IRT Model with Quadratic Effects
###########################################
#JAGS code for model with fixed evolution parameter
dyn_irt_model_quadratic <- 'model{
#loop through actors
for(i in 1:nactors){
#loop through votes
for(j in 1:nvotes){
Y[i,j] ~ dbern(mu[i,j])
logit(mu[i,j]) <- alpha[j] + beta[j] * theta[actor[i], period[j]] + beta2[j] * theta[actor[i], period[j]] * theta[actor[i], period[j]]
}
change[i] <- theta[i,2] - theta[i,1]
}
#set normal priors on betas and alphas
for(j in 1:nvotes){
beta[j] ~ dnorm(0, 1)
beta2[j] ~ dnorm(0, 1)
alpha[j] ~ dnorm(0, 1)
}
#set priors on thetas (fix space with two actors)
theta[1, 1] ~ dnorm(1, 1)
theta[2, 1] ~ dnorm(-1, 1)
for(t in 2:nperiods){
theta[1, t] ~ dnorm(theta[1, t-1], tau.evol)
theta[2, t] ~ dnorm(theta[2, t-1], tau.evol)
}
for(c in 3:nactors){
theta[c, 1] ~ dnorm(0, 1)
for(t in 2:nperiods){
theta[c, t] ~ dnorm(theta[c, t-1], tau.evol)
}
}
#set priors on evolution parameters for thetas
tau.evol <- 1
}
'
#Set parameters
jags.data.irt <- list(
Y = rolls.jags,
nactors = 630,
nvotes = 50,
nperiods = 2,
period = c(2, 2, 2, 2, 2, 2, 2, 2, rep(1, 42)),
actor = as.numeric(c(1:630))
)
#Fit JAGS model
model.jags.irt <- jags.model(file=textConnection(dyn_irt_model_quadratic), data = jags.data.irt, n.chains = 1, inits=list(.RNG.name="base::Wichmann-Hill",.RNG.seed=2355328))
posterior.coda.irt <- coda.samples(model.jags.irt, c("theta","beta","beta2","alpha", "tau.evol", "change"), n.iter=10000, thin = 10, progress.bar="text")
#Extract debate effects (betas)
sum.jags.irt=summary(posterior.coda.irt)[[1]][,1]
sum.jags.irt.beta=data.frame(means=sum.jags.irt[grepl(names(sum.jags.irt),pattern="beta")])
sum.jags.irt.beta$debate=as.numeric(gsub("\\D","",row.names(sum.jags.irt.beta)))
#Extract actor positions (thetas)
sum.jags.irt.theta=data.frame(means_period1_sq=sum.jags.irt[grepl(names(sum.jags.irt),pattern="theta")&
grepl(names(sum.jags.irt),pattern=",1")],
means_period2_sq=sum.jags.irt[grepl(names(sum.jags.irt),pattern="theta")&
grepl(names(sum.jags.irt),pattern=",2")],
change_sq=sum.jags.irt[grepl(names(sum.jags.irt),pattern="change")])
#Get 95% credible intervals
brexit.irt=HPDinterval(posterior.coda.irt[[1]],prob=.95)
#Save CIs in data frame
brexit.irt.thetas=data.frame(lower_period1_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "theta")&grepl(row.names(brexit.irt),pattern = ",1"),"lower"],
upper_period1_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "theta")&grepl(row.names(brexit.irt),pattern = ",1"),"upper"],
lower_period2_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "theta")&grepl(row.names(brexit.irt),pattern = ",2"),"lower"],
upper_period2_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "theta")&grepl(row.names(brexit.irt),pattern = ",2"),"upper"],
lower_change_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "change"),"lower"],
upper_change_sq=brexit.irt[grepl(row.names(brexit.irt),pattern = "change"),"upper"])
#Add MP names and actor numbers
brexit.irt.thetas$MP=rolls$mpname
brexit.irt.thetas$MPnumber=1:nrow(rolls)
sum.jags.irt.theta$MP=rolls$mpname
sum.jags.irt.theta$MPnumber=1:nrow(rolls)
#Join CIs and means
brexit.irt.thetas_sq=left_join(brexit.irt.thetas, sum.jags.irt.theta)
#Assess convergence of model with Geweke diagnostic
geweke_irt <- unlist(geweke.diag(posterior.coda.irt))
geweke_irt.dat=data.frame(gew=geweke_irt[!grepl("change",names(geweke_irt))&!grepl("frac",names(geweke_irt))],
names=names(geweke_irt[!grepl("change",names(geweke_irt))&!grepl("frac",names(geweke_irt))]))
geweke_irt.dat=geweke_irt.dat%>%filter(!is.na(gew))
#Assess how many parameters fall in 2/-2 range
prop.table(table(geweke_irt.dat$gew>2|geweke_irt.dat$gew< -2))
table(geweke_irt.dat$gew>2|geweke_irt.dat$gew< -2) #looks as expected - good convergence
#Plot Geweke diagnostics
geweke_irt_graph=ggplot(geweke_irt.dat) +
geom_histogram(binwidth=0.2,aes(x=gew,y=..density..), position="identity") +
geom_density(aes(x=gew,y=..density..))+
labs(x="Geweke Statistics",y="Density")+
theme_bw()+
theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.text=element_text(size=12,colour="black",family="Verdana"),
axis.title=element_text(size=12,family="Verdana"),
legend.position="none")+
scale_x_continuous(limits = c(-4,4))
geweke_irt_graph #looks like convergence
#Rename data
MP_data <- brexit.irt.thetas_sq
#Export data
write.csv(MP_data, "./Created auxiliary datasets/irt_models_MPs.csv")
